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Abstract 

We present a numerical study of Riemann's formula for the oscillating part of the density of the 
primes and their powers. The formula is comprised of an infinite series of oscillatory terms, one 
for each zero of the zeta function on the critical line and was derived by Riemann in his paper 
on primes assuming the Riemann hypothesis. We show that high resolution spectral lines can 
be generated by the truncated series at all powers of primes and demonstrate explicitly that the 
relative line intensities are correct. We then derive a Gaussian sum rule for Riemann's formula. 
This is used to analyze the numerical convergence of the truncated series. The connections to 
quantum chaos and semiclassical physics are discussed. 
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I. INTRODUCTION 



The Riemann zeta function is at the heart of number theory, and has also played a pivotal 
role in the study of quantum chaos Q| . There is a deep connection between the complex zeros 
of the zeta function and Random Matrix Theory j^. The zeros possess the same statistical 
properties as the energy eigenvalues of a dynamical Hamiltonian that is nonintegrable and 
whose dynamics are not time-reversal invariant. Unfortunately, this Hamiltonian is not 
known in terms of its dynamical variables. The main source of insight into this unknown 
quantum chaotic system comes from Gutzwiller's pioneering work [3], that expresses the 
oscillatory part of the quantum density of states as a sum over classical periodic orbits. 
(Such sums are now referred to as trace formulas.) It is well-known that the oscillatory part 
of the density of the Riemann zeros is given by a Gutzwiller-like sum, with one periodic term 
for every integer power of a prime number (A smoothed density of the Riemann zeros 
has also been studied in Ref. j^.) From this perspective, one can infer that a spectrum 
consisting of the Riemann zeros is generated by a Hamiltonian (albeit unknown) whose 
classical orbits have actions that are logarithms of primes and integer powers of primes. 

Conversely, one could ask whether it is also possible to calculate the prime number se- 
quence from a sum of oscillatory terms, with one term for every zero of the zeta function. 
Although less widely-known, such a series was actually given by Riemann himself 0]. Rie- 
mann derived an exact formula for the density of the primes (and their integer powers) that 
can be expressed as the sum of a smooth function and an infinite series of oscillatory terms 
involving the complex zeros of the zeta function. The smooth part has been thoroughly 
studied in the context of the prime number theorem whereas the oscillatory part has been 
largely ignored. Interestingly, it is the latter that contains the essential information about 
the location of the primes, as shown below. There is a vast literature on the distribution 
of the prime numbers. It is recognized that their distribution exhibits global regularity and 
local irregularity The nearest-neighbour spacings (NNS) of the primes is known to be 
Poisson-like j^, corresponding to an almost uncorrelated random distribution. This is very 
different in character from the Gaussian Unitary Ensemble (GUE) distribution of the Rie- 
mann zeros. Nevertheless, it is possible to generate the almost-uncorrelated sequence of all 
integer powers of primes from the interference of the highly-correlated Riemann zeros. 

As mentioned above, from the perspective of semiclassical periodic orbit theory, the 
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density of the Riemann zeros has the structure of a dynamical trace formula with periodic 
orbits. It is natural to ask whether Riemann's formula is a trace formula for the primes. 
Despite having the oscillatory terms, as discussed below, Riemann's formula is not a trace 
formula of dynamical origin. But, this does not preclude the existence of a trace formula 
for the primes. If one does exist, then this would support the notion that there exists a 
Hamiltonian system whose quantum spectrum is the primes. In any case, the exclusion 
of Riemann's formula as a trace formula suggests that there would be no correspondence 
between the classical dynamics and the Riemann zeros for this system. 

The purpose of this paper is to study the density of the primes from the perspective of 
numerical semiclassics. To our knowledge, Riemann's scries has not been studied numeri- 
cally. We first verify that Riemann's formula does produce spectral lines at the positions 
of the primes and their powers, even when the series is truncated. This is not completely 
unexpected since Riemann's series converges conditionally to the exact density which is a 
set of (5-function spikes. However, the 5 functions arise from the entire series. The truncated 
series is an approximation to the exact density. It does not yield spikes, but rather lines of 
various widths, heights and (unknown) shapes and it is not at all obvious that the relative 
line intensities of the truncated series are correct. Wc examine this problem both numeri- 
cally and analytically. We then provide a simple rule for estimating the value of the largest 
zero required to sufficiently resolve individual lines of a specific shape in some interval of 
interest and describe how to control the error incurred from a truncation of the series. 



II. RIEMANN'S FORMULA 



We start from the Euler product formula 

CW) = ll(^-P~T\ Re/5>1, (1) 
p 

where the product is over all primes p. It follows that X^^^i n exp(— Inp) = ln({f3). 
Dividing both sides by (3 and then taking the inverse Laplace transform of both sides with 
respect to E, we immediately obtain 

^(^) = E E - 1^^^") = ^ ^e/^^d/3, (2) 
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where a > 1. Riemann evaluated the RHS of Eq. Q to obtain N{E). Upon differentiation 

with respect to E and the subsequent substitution x = e^, we obtain the density p{x) of 

p"'/n along the real axis x as 

, , 1 1 cosfo; Inx) . . 

Pix) = — TTT 2 > — ^, (3 

Inx x x^ — 1 mx x^'"^mx 

where x > 1. This formula assumes the Riemann hypothesis, which states that the infinite 
number of complex zeros of the zeta function all lie on the critical line (3 = (1/2 ± ia), 
where a is real and positive. Note that explicit use of the symmetry of the complex zeros 
has been made to reduce the summation to cosine functions. A generalized version of the 



Riemann formula, where the zeros may lie anywhere in the critical strip is given in Ref. 
We shall denote the sum over the oscillatory terms on the RHS of Eq. (jH)) as p(x). Since 
Eq. © is exact, it is clear that the 5-function spikes of p(x) must be generated from the 
interference of the terms in p(x). From our experience in periodic orbit quantization jlol |. 
we know that a coarse-grained version of the exact density of states can be reproduced even 
from a truncated periodic orbit sum. Therefore, in the following section, we focus on the 
numerical analysis of the truncated series. 

Before presentin g th e results, however, we briefly review the pioneering numerical work 
of Riesel and Gohl The LHS of Eq. Q is a set of step functions, with unit steps at 
every prime p, one-half steps at p^, one-third steps at and so on, and may be obtained by 
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taking the contour integral of \n({(3)/l3 on the RHS of the equation. Riemann 
this function by /(x) and Edwards by J(x) The number of primes less than x, denoted 
by 7r(x), may be expressed in terms of /(x) as 

= f ^(")/'-'^"). (4) 

n=l 

where fi{n) is the Mobius function The modulating effect of the oscillatory terms due 
to the first twenty nine pairs of the complex Riemann zeros was numerically examined by 



Riesel and Gohl ll| in 1970. This early work already showed the approximate formation of 
the first few steps at the prime numbers, and modulations for some larger primes. Note that 
the series dH) requires a knowledge of the Mobius function, and is much more complicated 
than Eq. Riesel and Gohl actually replaced the sum over the Mobius functions by the 
Gram series involving factorials which are difficult to compute accurately for large integers. 
In this paper, we shall rather study formula Q for p(x) since it contains more information 



than the formula for the density of the primes (no powers) which can be obtained from the 
derivative of 7i{x). 



III. NUMERICS 



Numerically, we can only evaluate a finite number of terms from Riemann's infinite series. 
Although it would seem by inspection that all terms of the series are equally important and 
that there is no optimal ordering of the terms, Riemann states that the series is conditionally 
convergent and that it must be summed in order of increasing size of a. (For any series whose 
convergence is conditional, the order of summation must be specified since different orderings 
produce different results.) In this case, the "natural order" is the correct one. Riemann goes 
on to state that with this ordering the truncated series should give an approximation to the 
density of primes (and their powers), but that using a different ordering, the resulting finite 
series can approach arbitrary real values. We have verified this numerically and found that 
using finite sets of zeros chosen according to different rules yields incorrect results. Thus, 
for the numerical work that follows, we use the correct ordering. 



A. Line Intensities of the Truncated Series 



We first computed Eq. © using the first 10^ zeros and observed lines at the positions 
of the primes and their powers for x < 5000. However, for x > 2000, many lines cannot 
be fully resolved and the signal eventually dies out. This is due to truncation since only a 
small number of zeros have been included. (This will be discussed in more detail below.) 
Nevertheless, even this small number of zeros yields narrow lines at the lowest primes. In 
Fig. ^ we display the result for x G [1.5, 100]. Although one can clearly observe lines at the 
positions of the primes, the relative intensities cannot be determined by inspection since the 
lineshapes are not uniform (see Fig. |21). This is a common problem in spectral analysis, and 
is often resolved by imposing a more uniform lineshape through convolution of the signal 
with an appropriate smooth "response function" The response function is typically a 
peaked function that falls to zero in both directions from its maximum. Gaussian functions 
are positive-definite and decay rapidly. They are also convenient to use since their shape 
only depends on a single parameter (the variance) and therefore can be easily controlled. 
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FIG. 1: The result of computing Eq. © using the first 10'^ zeros for x G [1.5, 100]. The inset shows 
a closer view of three lines that appear for x £ [77, 84] . The three lineshapes are similar so that their 
relative heights are somewhat meaningful. However, the lineshapes vary considerably throughout 
the entire range so that heights cannot be immediately interpreted as relative intensities. 
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FIG. 2: A closer view of two non-adjacent lines in Fig. ^ The lineshapes clearly differ so that 
their relative heights are not meaningful. 

Thus, we next convolve this approximate density p{x) (i.e. smooth term and truncated 
series) with an unnormalized Gaussian of variance a: 

/oo 
p{x')G„{x — a;')dx', (5) 

where 

G„{x) = exp{-xy2a^). (6) 
6 
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The effect of the convolution is that rapidly oscillating features are washed out and smooth 
peaked features are smeared into the shape of the response function. If the lines were perfect 
S functions of height -D„, then from Eq. these would be replaced by DnG„{x — 
i.e. Gaussians of variance cr with height at a; = x„. Of course, the lines are not 5- 
function spikes so that the resultant lineshapes are not exactly Gaussian, but as long as the 
intrinsic linewidth is sufficiently small compared to the variance, the deviation from a perfect 
Gaussian is quite negligible. Therefore, the convolution produces a series of Gaussian lines, 
each of the same width. The key point is that the lineshapes are now essentially uniform so 
that the actual heights can be meaningfully compared and immediately interpreted as the 
relative intensities. It is important to keep in mind that since the response function has a 
maximum height of unity, the height of a line after convolution should be the area under that 
line before convolution. The reason for this is that although the lines of the original signal 
act like S functions with respect to the response function, they do have nonzero widths 
and so their effective 5-function "heights" Dn, are equal to the areas. In this sense, the 
convolution procedure is equivalent to directly integrating the area under each line of the 
signal. However, the convolution technique is much simpler and avoids errors that can arise 
from the long oscillatory tails of individual lines. Note that this procedure cannot resolve 
two adjacent lines when the spacing between them is smaller than a and thus o"max = 1/2. 

We have computed Eq. (jSJ for the range of interest in Fig. ^ This is shown in Fig. El 
using a = 0.05. We note here that the heights do not depend on the specific value of a 
due to the fact the Guassian is not normalized. As more terms are included in the sum, 
the natural linewidths decrease and the convolution becomes more accurate. It is then 
possible to produce high-resolution lines by using smaller variances. For example, using 
10^ zeros, we produced lines with a variance cr = 0.01. It would be useful to know how 
many primes can be resolved using a prescribed number of zeros. In the present scheme, 
one simply observes where the lines of the original signal develop a sufficiently large width. 
The important criterion here is that all linewidths should be at least smaller than the mean 
spacing between all powers of primes in the interval of interest. Of course, the width of 
any line is related to the number of terms used in truncating the series. Although this 
relationship can be determined, there is still the problem that all the lines have different 
shapes. Thus, we shall find it more useful to determine for what values of x the truncated 
formula can no longer produce lines of a specific uniform shape. 
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FIG. 3: The result of computing Eq. using the first 10^ zeros and a = 0.05. The range is the 
same as in Fig. ^ Note that hues with height less than unity occur at powers of primes and 
have height 1/n (for example, the line at 2^ = 64 has height 1/6 = 0.16). 

B. Sum Rule and Numerical Convergence 

The convergence of the series can be examined using a more controlled application of 
the convolution procedure described above. The general idea is to construct a series ( "sum 
rule") that absolutely converges to a "coarse-grained" version of the exact density. This 
density is obtained by replacing all spikes of the exact density by smooth peaked functions. 
One immediate advantage is improved convergence since it is easier to reproduce the well- 
defined smooth peaked functions of a coarse-grained density using a truncated sum rule than 
it is to reproduce spikes using the original truncated series. However, for our purposes, the 
more important reason for using a sum rule is to control the convergence of the series. This 
will become evident after the sum rule is given. The sum rule itself is obtained from a direct 
convolution of the original series with a "smoothing function" , that is, some smooth function 
whose Fourier coefficients rapidly decrease. (Since the original series consists of cosines, the 
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resulting integral is essentially a cosine transform of the smoothing function.) 

The above discussion is quite general. We now connect this idea with the numerical 
calculations described above. Assume the coarse-grained density consists of a set of Gaussian 
functions of variance a centered at each prime (or prime power) with heights equal to unity 
(or the reciprocal power). We want to construct a series as described above that converges 
to this density, that is, we want to find a Gaussian sum rule for Riemann's series. To do 
this, we convolve Riemann's series term-by-term with a Gaussian smoothing function. For 
the following calculation, we will define Sa{x) = alnx and A{x) = — 2/y^lnx. Then, we 
can write Riemann's formula as p{x) = A{x) Re{exp[iS'a(x)]}. The Gaussian sum rule 
is 



p^{x) = ~p{x) * G„{x) = / p(x')e-("-"') /^a^dx' 

J —oo 

= J]Re ( r A(a;')e^^"(^')e-(^-^')'/2-^dx' |> , a; > 1. (7) 



For a < o"max = 1/2, the Gaussian rapidly decays to zero. This implies that the main 
contribution to the integral comes from a small interval centered about x' = x. Elsewhere, 
the integrand is practically zero. Thus, we make two approximations to proceed further. 
First, the amplitude function A{x') changes very slowly and on the small interval of interest 
A{x') ~ A{x). Secondly, the phase function Sa{x') can be replaced by its Taylor series 
expansion about x' = x: Sa{x') = Sa{x) + S'^{x) {x' — x) + . . .. If we retain the leading-order 
term only. 



p^{x) ^ A(a;)^Re| e*[5.(-)+5;W{x'-x)]g-{x-x')V2<x2d^/ 

= A{x)^Re |e^[Sa{x)-xS;(x)]g-xV2a2 ^-lx'^-i2x+2ia^S'M)W]/^^' (^^'^ 

= \/2^aA(x)5^e-'^'^^'(")/2j^e{e*^"(")}, (8) 

a 

where we have used the standard result for the Gaussian integral Finally, the Gaussian 
sum rule for Riemann's series is 

p^{x) = -^^^y e-'^'"'/2.2cos(alnx). (9) 
wx m X ^-^ 

This sum rule explicitly shows the effect of convolution on the series; each term is modulated 
by an exponential factor. This factor essentially controls the convergence of the series for 



all values of x. Although the orginal series is only conditionally convergent, as long as the 
correct ordering is used, this sum rule is also absolutely convergent. As stated above, we seek 
an approximate relation between the maximum zero included in the sum and the maximum 
prime that can be resolved. One way to determine this is as follows. First, specify the value 
of the largest zero, Omax and include all zeros a < ttmax- Then, there exists a set of values 
X < Xmax for which the exponential factor falls below some threshold parameter e. This 
condition immediately gives the simple relation 



where < e < M. For a > Omax and x < Xmax? all terms are exponentially smaller than 
e and are thus numerically insignificant. The choice of the parameter e depends on the 
desired precision of a resolved line. An upper bound M for the par ameter is the value of the 
exponential factor (e~^/^) at its inflection point xi = {l/\/3)aa [l5|. This implies Xmax < xj. 
The lower bound can be as small as machine zero (for example 10~^^). However, there is 
no reason for such an extreme choice since we are mostly interested in determining where 
numerical errors become significant (i.e. where lines are no longer visibly resolved and the 
intensities are erroneous by more than 1%). Of course, higher precision can be imposed at 
the cost of resolving fewer primes. But, since the improved precision will not be apparent 
in the graph of Pa{x), there is no compelling reason to choose exceedingly small values. For 
our purposes, a convenient choice is e = e~^/^. 

We now provide a few examples to illustrate the utility of relation (fTIUl . As a first example, 
we take Omax = 9878 which is the lO^th zero. Using the above formula (with a = 0.1) yields 
a^max = 373. In Fig. HI we evaluate Eq. Q using the first 10"^ zeros (and include the smooth 
term). One can clearly see significant errors for x > 400. As a second example, we take 
ttmax = 74921 (the lO^th zero). The formula then gives Xmax = 2832. In Fig. we truncate 
Eq. (jni) at this value of amax and observe significant errors occur for x > 2900. 

An additional benefit of the sum rule is that it gives us an immediate measure of the 
error incurred from truncation. The largest errors are in the vicinity of Xmax where there are 
contributions 0{e) that have been excluded. For all other values of x < Xmax, the excluded 
terms are exponentially smaller. Of course, we have complete control of this error through 
our freedom in specifying e. In the case of the original truncated series, it is not immediately 
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FIG. 4: The sum rule Q using the first 10^ zeros and a = 0.1. The lower window shows a closer 
view of the lines in the interval 0.99 <Z pcr{x^ <Z 1. Formula (jlOj) indicates that for x ^ ^^max ~ 

373, 

all lines should be resolved and intensities should have errors less than a"^ = 0.01. 

obvious what the errors are, but they can be determined through more elaborate analysis. 



IV. DISCUSSION AND CONCLUSION 



By writing (^(7 + it) = + it) \ exp(— i^^(t)), we see that all the information about the 
zeros along the t— axis is contained in the phase 6^{t). This has to jump by it to accomodate 
the sign change in ( at every zero, and it can be shown that the oscillating part of the 
density of the zeros on the critical line is proportional to the derivative of the imaginry part 
of InC(t) with respect to t On the other hand, we see from Eq. Q that the appropriate 
contour integral over \n({f3) also yields p{x) relating to the primes. Thus, the phase of the 
zeta function, as defined above, connects the Riemann zeros to the primes. 

As mentioned above, if the series is truncated, the signal gradually dies out as x increases. 
This can be understood by noting that due to the logarithmic dependence, each term pro- 



11 



0.8 
„ 0.6 

X 

=^ 0.4 
0.2 


2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 
X 



1 

0.999 
0.998 
0.997 
^ 0.996 
0.995 
0.994 
0.993 
0.992 
0.991 
0.99 

2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 
X 



FIG. 5: The sum rule Q using the first 10^ zeros and a = 0.1. The lower window shows a closer 
view of the lines in the interval 0.99 <Z p^ji^x^ <z 1. Formula (jlOjl indicates that for x ^ ^max ~ 

2832, 

all lines should be resolved and intensities should have errors less than o"^ = 0.01. 



duces an oscillation whose period continually increases while its amplitude decays. Clearly, 
more high frequency (large a) terms are required for sufficient constructive interference. 
This explains the fact that lines at small values of x are resolved more quickly than at larger 
values. Although the higher frequency terms are responsible for short-range oscillations 
and one could imagine exclusive use of those terms rather than lower frequency terms, the 
difficulty is the conditional convergence of the series and the fact that all of the terms are 
equally important. Unfortunately, this implies that Riemann's formula is impractical for 
resolving lines at large primes. This is also consistent with Eq. (fTUI) . If one is interested 
in using Riemann's series to find (new) large primes, for example, on the order of 10^^°'°°'^, 
then one requires an accurate knowledge of roughly the same number of zeros. 

We emphasize that Riemann's formula, is only correct if the Riemann hypothesis is true. 
Otherwise, if a pair of zeros occur at f3± = 'y ± ia, the factor x^^^ in the denominator of 
the oscillating term of Eq. Q should be replaced by x^^~'^^ P]. An interesting numerical 
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experiment is to move the zeros off the critical hne, that is, to arbitrarily change their real 
parts. We find this still produces lines at the primes and their powers, but the relative 
intensities are incorrect. This is interesting since it demonstrates that the location of the 
primes depends only on the imaginary part of the zero. The real part only affects the 
intensities which are anyway not evident from a direct evaluation of the series. This provides 
another motivation for the numerical and analytical procedures described in this paper. 

It is natural to compare the oscillating part of the density p{x) with the semiclassical 
trace formula ^, [3| of a dynamical system. One could identify a as an orbit label, one 
for each zero of the zeta function and x as the single-particle energy variable. Then, p{x) 
in Eq. may be interpreted as the density of states as a function of energy with the first 
term on the RHS corresponding to the smooth Thomas- Fermi (TF) contribution In the 
oscillating part, the argument a In x of the cosine term should then correspond to the action 
Sa{x) of the orbit a. Note however that there are no implicit repetition indices in Eq. Q 
thereby implying that even if one gives a dynamical interpretation to p{x), the orbits are 
not periodic. This is in direct contrast to the trace formula for the Riemann zeros, in which 
the orbits are periodic with primitive period \np for each prime j^. Of course, the most 
striking feature is that the amplitude has no a dependence. Even oscillatory contributions 
to the density of states from non-periodic trajectories usually have amplitudes that depend 
on the orbit In the event that there is a fortuitous cancellation of the index a, it is 
unlikely that the energy dependence in the denominator of the oscillating term as well as 
the TF term can then be generated consistently by the same Hamiltonian. Consequently, 
Riemann's formula is not a trace formula of dynamical origin. 

With regard to spectral statistics, it is well known that nearest-neighbour spacings (NNS) 



19| of the Riemann zeros obey the GUE distribution of Random Matrix 
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'heory, character- 



2l\ . The same zeros 



istic of a chaotic quantum system without time-reversal symmetry 
also generate the spectrum of the primes and their powers through Riemann's formula Q. 
As mentioned earlier, the NNS distribution of the primes is Poisson-like , with some level 
repulsion, which, if at all of dynamic origin, hints only to near-integrability Q]. Thus, it is 
quite remarkable that the highly-correlated sequence of the zeros can (through Riemann's 
formula) interfere to produce the almost-uncorrelated sequence of the primes. 

In conclusion, we have demonstrated that the spectrum of the primes and their integer 
powers can be accurately generated from a sum of periodic terms, each term involving a 
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zero of the zeta function. This is in the spirit of semiclassical periodic orbit theory, where 
the individual levels of a quantum spectrum may be resolved from a sum of oscillatory 
terms, each arising from periodic orbits. Despite the accuracy of the generated spectrum, 
Riemann's formula is not a trace formula. However, this does not imply that there is no 
such formula, and it would still be interesting to understand the spectrum of the primes in 
terms of periodic orbits. This could provide insight into the structure of a possible trace 
formula for the primes. If this formula could be found, the remaining challenge would be to 
obtain the corresponding Hamiltonian. 
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